! ---
! Copyright (C) 1996-2016	The SIESTA group
!  This file is distributed under the terms of the
!  GNU General Public License: see COPYING in the top directory
!  or http://www.gnu.org/copyleft/gpl.txt .
! See Docs/Contributors.txt for a list of contributors.
! ---
module m_digest_nnkp

  use precision,         only: dp

  implicit none

  integer              :: iun_nnkp    ! Integer that points to the .nnkp file
                                      !   used in read_nnkp and  scan_file_to

CONTAINS
! subroutine read_nnkp
! subroutine scan_file_to
! function   vectorproduct
! subroutine chosing_b_vectors
! subroutine set_excluded_bands
! function   getdelkmatgenhandle
! subroutine number_bands_wannier


subroutine read_nnkp( seedname, latvec, reclatvec, numkpoints,          &
                      kpointsfrac, nncount, nnlist, nnfolding,          &
                      numproj, projections, numexcluded, excludedbands, &
                      w90_write_amn )
!-----------------------------------------------------------------------
! This subroutine reads the .nnkp file generated by the Wannier90 code
! run in preprocessor mode, and populates some variables defined in the
! siesta2wannier module.
! The meaning of all these variables are properly described
! in the file siesta2wannier90.F90 file
!-----------------------------------------------------------------------
  use sys,               only: die               ! Termination routine
  use siesta_geom,       only: ucell             ! Unit cell lattice vectors
  use units,             only: Ang               ! Conversion factor
  use alloc,             only: re_alloc          ! Reallocation routines
  use m_mpi_utils,       only: broadcast         ! Broadcasting routines
  use files,             only: label_length      ! Number of characters allowed
                                                 !   in the name of a file
  use parallel,          only: IOnode            ! Input/output node
  use parallel,          only: Node              ! Working node
  use trialorbitalclass                          ! Derived class to define the
                                                 !   trial wave functions

  implicit none

  character(label_length+3),       intent(in)  :: seedname ! Name of input file
  integer,                         intent(out) :: numkpoints
  integer,                         intent(out) :: nncount
  integer,                         intent(out) :: numproj
  integer,                         intent(out) :: numexcluded
  integer,  pointer       :: nnlist(:,:)
  integer,  pointer       :: nnfolding(:,:,:)
  integer,  pointer       :: excludedbands(:)
  real(dp),                        intent(out) :: latvec(3,3)
  real(dp),                        intent(out) :: reclatvec(3,3)
  real(dp), pointer       :: kpointsfrac(:,:)
  type(trialorbital), allocatable, intent(out) :: projections(:)
  logical,                         intent(out) :: w90_write_amn

  external :: io_assign, io_close

!
! Internal variables
!
  integer              :: i, j, ik, iw! Loop counters
  integer              :: ikp, inn    ! Loop counters
  integer              :: iikp        ! Dummy variable
  integer              :: iu          ! Integer number to refer to the nnkp file
  logical              :: have_nnkp   ! Yes if the .nnkp is found
                                      ! No if the .nnkp file is not found
  real(dp)             :: rcell(3,3)  ! Reciprocal lattice vectors
! Variables related to the projection functions.
  real(dp)             :: center_w(3) ! Projection function center, 
                                      !   as read from the nnkp file
                                      !   (relative to the direct lattice
                                      !   vectors)
                                      !   They are transformed latter to 
                                      !   cartesian coordinates.
  integer              :: l_w         ! Angular momentum of the projection func.
  integer              :: mr_w        ! z-projection quantum number of the
                                      !   projection function
  integer              :: r_w         ! Radial quantum number of the 
                                      !   projection function
  real(dp)             :: zaxis(3)    ! Defines the axis from which the polar
                                      !   angle theta in spherical 
                                      !   polar coordinates is measured.
                                      !   Default: (0.0 0.0 1.0)
  real(dp)             :: xaxis(3)    ! Defines the axis from which the
                                      !   azimuthal angle phi in spherical 
                                      !   polar coordinates is measured.
                                      !   Must be orthogonal to z-axis.
  real(dp)             :: zovera      ! z/a, diffusivity, spread.
                                      !   Units in nnkp file: Ang^-1.
                                      !   This is transformed latter to Bohr^-1
  real(dp)             :: xnorm       ! Norm of the x-vector 
  real(dp)             :: znorm       ! Norm of the z-vector 
  real(dp)             :: coseno      ! Cosine between the x and z vectors

  real(dp), parameter :: eps4  = 1.0e-4_dp
  real(dp), parameter :: eps6  = 1.0e-6_dp
                                       

! Initialize some variables
  nullify( kpointsfrac   )
  nullify( nnlist        )
  nullify( nnfolding     )
  nullify( excludedbands )

! Check the existence of the .nnkp file and open it if it exists.
  if (IOnode) then  

     inquire(file=trim(seedname)//".nnkp",exist=have_nnkp)
     if( .not. have_nnkp ) then
       write(6,'(/,a)')  &
 &       'read_nnkp: Could not find the file '//trim(seedname)//'.nnkp'
       call die()
     endif

     call io_assign( iu )
     iun_nnkp = iu
     open( iun_nnkp, file=trim(seedname)//".nnkp",form='formatted')

  endif ! IOnode

  if (IOnode) then  
!   check the information from *.nnkp with the nscf_save data
    write(6,'(/,a)')  &
 &    'read_nnkp: Checking info from the ' // trim(seedname) // '.nnkp file'

!
!   Real lattice vectors 
!
    call scan_file_to('real_lattice')
       write(6,'(a)')'read_nnkp: Reading data about real lattice'
       do j = 1, 3
         read(iun_nnkp,*) (latvec(i,j),i=1,3)
!        Rescale from Ang to Bohr, that are the units internally used in Siesta
         do i = 1, 3
           latvec(i,j) = latvec(i,j)*Ang
         enddo
       enddo
!      Compare the unit cell read from .nnkp with the unit cell used in Siesta
       do j = 1, 3
         do i = 1, 3
           if(abs(latvec(i,j)-ucell(i,j))>eps4) then
             write(6,*)  'read_nnkp: Something wrong with the real lattice! '
             write(6,*)  ' latvec(i,j) =',latvec(i,j),' ucell(i,j)=',ucell(i,j)
             call die()
           endif
         enddo
       enddo
       write(6,'(a)')'read_nnkp:      - Real lattice is ok'

!
!   Reciprocal lattice vectors 
!
    call scan_file_to('recip_lattice')
       write(6,'(a)')'read_nnkp: Reading data about reciprocal lattice'
       do j = 1, 3
         read(iun_nnkp,*) (reclatvec(i,j),i=1,3)
!        Rescale from Ang^-1 to Bohr^-1, that are the units internally 
!        used in Siesta
         do i = 1, 3
           reclatvec(i,j) = reclatvec(i,j)/Ang
         enddo
       enddo 
!      Compute the reciprocal lattice from the unit cell used in Siesta
       call reclat(ucell, rcell, 1)
!      Compare the reciprocal lattice read from .nnkp with the one used inSiesta
       do j = 1, 3
         do i = 1, 3
           if(abs(reclatvec(i,j)-rcell(i,j))>eps4) then
             write(6,*)'read_nnkp: Something wrong with the reciprocal lattice!'
             write(6,*)' reclatvec(i,j)=',reclatvec(i,j), &
 &                     ' rcell(i,j)=',rcell(i,j)
             call die()
           endif
         enddo 
       enddo 
       write(6,'(a)')'read_nnkp:      - Reciprocal lattice is ok'
  endif ! IOnode
! Broadcast the lattice vectors in real space and the reciprocal lattice vectors
  call broadcast( latvec     )
  call broadcast( reclatvec  )

!
! k-points
!

  if (IOnode) then  
    call scan_file_to('kpoints')
       write(6,'(a)')'read_nnkp: Reading data about k-points'
       read(iun_nnkp,*) numkpoints
  endif

! Broadcast information regarding the number of k-points 
! and allocate in all nodes the corresponding arrays containing information 
! about the k-points.
  call broadcast( numkpoints )
  call re_alloc(  kpointsfrac, 1, 3, 1, numkpoints,  &
 &                name='kpointsfrac', routine='read_nnkp')

  if (IOnode) then
     do ik = 1, numkpoints
       read(iun_nnkp,*) kpointsfrac(1:3,ik)
     enddo
  endif 

! Broadcast the information regarding the k-points
  call broadcast( kpointsfrac )

!
! Projections
!
  numproj = 0
  if (IOnode) then  ! Read nnkp file on ionode only
    write(6,'(a)')'read_nnkp: Reading data about projection centers '
    call scan_file_to('projections')
    read(iun_nnkp,*) numproj
  endif ! IOnode
  call broadcast( numproj )

  if( numproj .ne. 0 ) then
    if (IOnode) then   ! read from ionode only
      allocate(projections(numproj))
      do iw = 1, numproj
        read(iun_nnkp,*) (center_w(i), i=1,3), &
 &                       l_w,                  & 
 &                       mr_w,                 &
 &                       r_w
        read(iun_nnkp,*) (zaxis(i), i = 1, 3), &
 &                       (xaxis(i), i = 1, 3), & 
 &                        zovera

!       Check that the xaxis and the yaxis are orthogonal to each other
        xnorm = sqrt(xaxis(1)*xaxis(1) + xaxis(2)*xaxis(2) + xaxis(3)*xaxis(3))
        if (xnorm < eps6) call die('read_nnkp: |xaxis| < eps ')
        znorm = sqrt(zaxis(1)*zaxis(1) + zaxis(2)*zaxis(2) + zaxis(3)*zaxis(3))
        if (znorm < eps6) call die('read_nnkp: |zaxis| < eps ')

        coseno = (xaxis(1)*zaxis(1) + xaxis(2)*zaxis(2) + xaxis(3)*zaxis(3)) / &
 &               xnorm/znorm
        if (abs(coseno) > eps6) &
          call die('read_nnkp: xaxis and zaxis are not orthogonal')
        if (zovera < eps6) &
          call die('read_nnkp: zovera value must be positive')

!       Now convert "center" from ScaledByLatticeVectors to Bohrs
        do i=1,3
          projections(iw)%center(i) = 0.0_dp
          do j=1,3
            projections(iw)%center(i) = center_w(j)*latvec(i,j) + &
 &                                      projections(iw)%center(i)
          enddo
        enddo

        projections(iw)%zaxis  = zaxis
        projections(iw)%xaxis  = xaxis
!       Compute the y-axis
        projections(iw)%yaxis = vectorproduct(zaxis,xaxis)
        projections(iw)%yaxis = projections(iw)%yaxis /     &
 &         sqrt(dot_product(projections(iw)%yaxis,projections(iw)%yaxis))

        projections(iw)%zovera = zovera*Ang ! To Bohr^-1
        projections(iw)%r      = r_w
        projections(iw)%l      = l_w
        projections(iw)%mr     = mr_w

!       Select the maximum angular momentum
        select case(l_w)
          case(0:3)
            projections(iw)%lmax = l_w
          case(-3:-1)
            projections(iw)%lmax = 1 ! spN hybrids
          case(-5:-4)
            projections(iw)%lmax = 2 ! spdN hybrids
          case default
            call die("Invalid l in read_nkp")
        end select

!       Further checks
        if ( mr_w .lt. 1 ) &
 &        call die("Invalid mr in read_nnkp")
        if ( r_w .gt. 3 .or. r_w .lt. 1 ) then
          call die("Invalid r_w in read_nnkp")
        elseif ( l_w .ge. 0 .and. mr_w .gt. 2*l_w+1 ) then
          call die("Invalid mr_w in read_nnkp")
        elseif ( l_w .lt. 0 .and. mr_w .gt. 1-l_w ) then
          call die("Invalid mr_w in read_nnkp")
        endif

!       Cut-off initialization (in Bohr)
        projections(iw)%rcut = cutoffs(r_w)/projections(iw)%zovera 

      enddo
    endif ! IOnode

!   Broadcast the information about the projection functions
    call broadcast_projections( )

  endif ! numproj ne 0
 
!! For debugging
!  if (IOnode) then   ! read from ionode only
!    write(6,*) 'Projections:'
!    do iw = 1, numproj 
!       write(6,'(3f12.6,3i3,f12.6)') &
! &       projections(iw)%center(1:3),     &
! &       projections(iw)%l,               &
! &       projections(iw)%mr,              &
! &       projections(iw)%r,               &
! &       projections(iw)%zovera
!    enddo
!  endif ! IOnode
!! End debugging

  if ( numproj .eq. 0 ) then
    if( IOnode ) then
      write(6,'(/,a)')'read_nnkp: WARNING: No projections specified in ' &
 &      //trim(seedname)//".nnkp"
      write(6,'(a,/)')'read_nnkp: Calculation of Amn omitted'
      w90_write_amn = .false.
    endif
  endif

!
! k-points neighbours
!
  if (IOnode) then   ! read from ionode only
    write(6,'(a)')'read_nnkp: Reading data about k-point neighbours '
    call scan_file_to('nnkpts')
    read(iun_nnkp,*) nncount
  endif

! Broadcast information regarding the number of k-points neighbours
! and allocate in all nodes the corresponding arrays containing information 
! about k-point neighbours 
  call broadcast( nncount )
  call re_alloc( nnlist, 1, numkpoints, 1, nncount,           & 
 &               name = "nnlist", routine = "read_nnkp" )
  call re_alloc( nnfolding, 1, 3, 1, numkpoints, 1, nncount,  &
 &               name = "nnfolding", routine = "read_nnkp" )

  if (IOnode) then
    do ikp = 1, numkpoints
       do inn = 1, nncount
          read(iun_nnkp,*) iikp,                                &
 &                         nnlist(ikp,inn),                     &
 &                         ( nnfolding(j,ikp,inn), j = 1, 3 )
       enddo
    enddo
  endif

! Broadcast information regarding the list of k-points neighbours
  call broadcast( nnlist    )
  call broadcast( nnfolding )

!
! excluded bands
!
  if (IOnode) then     ! read from ionode only
    write(6,'(a)')'read_nnkp: Reading data about excluded bands '
    call scan_file_to('exclude_bands')
    read (iun_nnkp,*) numexcluded
  endif

! Broadcast information regarding the number of excluded bands
! and allocate in all nodes the corresponding arrays containing 
! information about excluded bands
  call broadcast( numexcluded )
  if( numexcluded .gt. 0 ) then
    call re_alloc( excludedbands, 1, numexcluded,     &
 &                 name = 'excludedbands', routine = "read_nnkp" )
  endif

  if (IOnode) then
    if( numexcluded .gt. 0 ) then
      do i = 1, numexcluded
        read(iun_nnkp,*) excludedbands(i)
      enddo
    endif
  endif

!! For debugging
!  if (IOnode) then
!    do i = 1, numexcluded
!      write(6,*) excludedbands(i)
!    enddo
!  endif
!! End debugging

! Broadcast information regarding excluded bands
  if( numexcluded .gt. 0 ) call broadcast( excludedbands )

  if (IOnode) call io_close(iun_nnkp)   ! ionode only

  return

end subroutine read_nnkp

!
!-----------------------------------------------------------------------
subroutine scan_file_to (keyword)
!
! This subroutine is used to digest the .nnkp file.
! It is borrowed from the QuantumEspresso suite.
!
!
  implicit none
  character(len=*)  :: keyword
  character(len=80) :: line1, line2
!
! by uncommenting the following line the file scan restarts every time
! from the beginning thus making the reading independent on the order
! of data-blocks
!   rewind (iun_nnkp)
!
10 continue
   read(iun_nnkp,*,end=20) line1, line2
   if(line1/='begin')  goto 10
   if(line2/=keyword)  goto 10
   return
20 write (6,*) keyword," data-block missing "
   call die("--")
end subroutine scan_file_to
!-----------------------------------------------------------------------

!
!-----------------------------------------------------------------------
function vectorproduct(a,b)
!
! This function computes the vector product of two vectors a and b
!
  real(dp),intent(in),dimension(3) :: a,b
  real(dp),dimension(3)            :: vectorproduct
  vectorproduct(1) = a(2) * b(3) - a(3) * b(2)
  vectorproduct(2) = a(3) * b(1) - a(1) * b(3)
  vectorproduct(3) = a(1) * b(2) - a(2) * b(1)
end function vectorproduct
!-----------------------------------------------------------------------

!
!-----------------------------------------------------------------------
subroutine chosing_b_vectors( kpointsfrac, nncount, nnlist, nnfolding, &
                              bvectorsfrac )
!
! In this subroutine we compute the b vectors that connect each
! mesh k-point to its nearest neighbours.
!
! ------------------------------ INPUT -----------------------------------------
! real(dp) kpointsfrac(:,:)  : List of k points where the overlap matrix
!                                of the periodic part of the wave function
!                                with its neighbour will be computed.
!                                Given in reduced units.
! integer nncount            : The number of nearest neighbours belonging to 
!                                each k-point of the Monkhorst-Pack mesh
! integer nnlist(:,:)        : nnlist(ikp,inn) is the index of the 
!                                inn-neighbour of ikp-point
!                                in the Monkhorst-Pack grid folded to the 
!                                first Brillouin zone
! integer nnfolding(:,:,:)   : nnFolding(i,ikp,inn) is the i-component 
!                                of the reciprocal lattice vector 
!                                (in reduced units) that brings
!                                the inn-neighbour specified in nnlist
!                                (which is in the first BZ) to the
!                                actual \vec{k} + \vec{b} that we need.
!                                In reciprocal lattice units.
! ------------------------------ OUTPUT ----------------------------------------
! real(dp) bvectorsfrac(:,:) : Vectors that connect each mesh k-point 
!                              to its nearest neighbours.
! ------------------------------ BEHAVIOUR -------------------------------------
! Set some kpointfrac, in this case the first in the list, and calculate
! b's by subtracting it from it's neighbors.
! They will be used as a handle to the "delkmat" matrix element of
! a plane wave with wave vector b.
! ------------------------------ UNITS -----------------------------------------
! bvectorsfrac are given in reduced units.
! ------------------------------------------------------------------------------

! Used module procedures
  use alloc,            only: re_alloc ! Re-allocation routine

  implicit none

  integer,  intent(in)           :: nncount
  integer,  pointer  :: nnlist(:,:)
  integer,  pointer  :: nnfolding(:,:,:)
  real(dp), pointer  :: kpointsfrac(:,:)
  real(dp), pointer  :: bvectorsfrac(:,:)

!
! Internal variables
!
  integer :: i, inn

  nullify( bvectorsfrac )
  call re_alloc( bvectorsfrac, 1, 3, 1, nncount,    &
                 name="bvectorsfrac", routine = "chosing_b_vectors")
! Set some kvector, in this case the first in the list and calculate
! b's by subtracting it from it's neighbors
! They will be used as a handle to the "delkmat" matrix element of
! a plane wave with wave vector b.
!
! Loop on the three spatial directions
  do i = 1, 3
!   Loop on neighbour k-points
    do inn = 1, nncount
      bvectorsfrac(i,inn) =                  &
        kpointsfrac(i,nnlist(1,inn)) + nnfolding(i,1,inn) - kpointsfrac(i,1)
    enddo
  enddo

end subroutine chosing_b_vectors


!
!-----------------------------------------------------------------------
integer function getdelkmatgenhandle( vec, nncount,  bvectorsfrac )
!
! This function accepts a k vector connecting two neighboring
! k-points and returns the position of vec in the array bVectorsFrac.
! It is used to assign matrix elements of a planewave to a k-vector.
! Vec must be in fractional coordinates.

  implicit none
  real(dp), intent(in) :: vec(3)
  integer,  intent(in) :: nncount                 ! Number of nearest k-point
                                                  !   neighbours  
  real(dp), intent(in) :: bvectorsfrac(3,nncount) ! Vectors that connect each 
                                                  !   mesh k-point to its
                                                  !   nearest neighbours.

! Internal variables
  integer              :: inn
  real(dp)             :: normSq,minNormSq

  minNormSq = 1.0_dp
  do inn = 1, nncount
    normSq = (bvectorsfrac(1,inn) - vec(1))**2 +                    &
 &           (bvectorsfrac(2,inn) - vec(2))**2 +                    &
 &           (bvectorsfrac(3,inn) - vec(3))**2
    if ( normsq .lt. minnormsq) then
      getdelkmatgenhandle = inn
      minnormsq = normsq
    endif
  enddo
end function getdelkmatgenhandle

!
!-----------------------------------------------------------------------
subroutine number_bands_wannier( numbandswan )

  use siesta_options,  only: hasnobup   ! Is the number of bands with spin up 
                                        !  for wannierization defined
  use siesta_options,  only: hasnobdown ! Is the number of bands with spin down
                                        !  for wannierization defined
  use siesta_options,  only: hasnob     ! Is the number of bands 
                                        !  for wannierization defined
                                        !  (for non spin-polarized calculations)
  use siesta_options,  only: nobup      ! Number of bands with spin up 
                                        !  for wannierization
  use siesta_options,  only: nobdown    ! Number of bands with spin down 
                                        !  for wannierization
  use siesta_options,  only: nob        ! Number of bands for wannierization
                                        !  (for non spin-polarized calculations)
  use m_spin,          only: nspin      ! Number of spin components
  use sparse_matrices, only: maxnh      ! Maximum number of orbitals
                                        !   interacting
  use sparse_matrices, only: numh       ! Number of nonzero elements of each
                                        !   row of the hamiltonian matrix 
  use sparse_matrices, only: listhptr   ! Pointer to start of each row of the
                                        !   hamiltonian matrix
  use sparse_matrices, only: Dscf       ! Density matrix in sparse form
  use sparse_matrices, only: S          ! Overlap matrix in sparse form
  use atomlist,        only: no_l       ! Number of orbitals in local node
                                        ! NOTE: When running in parallel,
                                        !   this is core dependent
                                        !   Sum_{cores} no_l = no_u
  use atomlist,        only: no_u       ! Number of orbitals in unit cell

  use m_noccbands,     only: number_of_occupied_bands
  use m_noccbands,     only: noccupied  ! Number of occupied bands for a given
                                        !   spin direction
  use parallel,        only: IOnode

#ifdef MPI
  use m_mpi_utils,     only: globalize_sum
#endif

  implicit none

! integer numbandswan(2) ! Number of bands for wannierization
  integer, intent(out), dimension(2) :: numbandswan 

! Internal variables
  integer :: ispin     ! Counter for the loops on the spin variables
  integer :: io        ! Counter for loops on the atomic orbitals
  integer :: j         ! Counter for loops on neighbours
  integer :: ind       ! Index of the neighbour orbital

  real(dp):: qspin(2)  ! Total population of spin up and down    

#ifdef MPI
  real(dp):: qtmp(2)   ! Temporary for globalized spin charge
#endif

! Initialize variables
  numbandswan(:) = 0

! Calculate the default number of bands considered for wannierization
! numbandswan = number of occupied bands
! This makes sense in molecules and insulators.
  do ispin = 1, nspin
    qspin(ispin) = 0.0_dp
    do io = 1, no_l
      do j = 1, numh(io)
        ind = listhptr(io) + j
        qspin(ispin) = qspin(ispin) + Dscf(ind,ispin) * S(ind)
      enddo
    enddo
  enddo

#ifdef MPI
! Global reduction of spin components
  call globalize_sum(qspin(1:nspin),qtmp(1:nspin))
  qspin(1:nspin) = qtmp(1:nspin)
#endif

!! For debugging
!   if( IOnode ) then
!     write(6,'(/,a,2f12.5)')  &
! &       'number_bands_wannier: qspin', qspin
!   endif
!! End debugging

  call number_of_occupied_bands( qspin )

! Setting up the default number of bands for wannierization to the 
! number of occupied bands
  do ispin = 1, nspin
    numbandswan(ispin) = noccupied(ispin)
  enddo

! Check the consistency with the required bands specified in the input file
  if ( nspin .eq. 1 ) then   ! Unpolarized calculation
!   If a number of bands is explicitly included in the input,
!   then use it overwriting the default value
    if ( hasnob ) then
      numbandswan(1) = nob
!     Since we are in a non-spin polarized calculation
!     ignore any possible value of the number of the number of bands for
!     wannierization introduced in the input
      if ( hasnobup ) then     
        if( IOnode )      &
 &        write(6,'(a)')  &
 &          'number_bands_wannier: ignoring Siesta2Wannier90.NumberOfBandsUp'
      endif
      if ( hasnobdown ) then
        if( IOnode )      &
 &        write(6,'(a)')  &
 &          'number_bands_wannier: ignoring Siesta2Wannier90.NumberOfBandsDown'
      endif
 
    elseif ( hasnobup .and. hasnobdown ) then
      numbandswan(1) = nobup
      numbandswan(2) = nobdown
      if( numbandswan(1) .ne. numbandswan(2) ) then
        if( IOnode )      & 
 &        write(6,'(a)')  &
 &          'number_bands_wannier: inconsistent output band specification'
        call die("Confusion in fdf")
      endif
    elseif ( hasnobup ) then
      numbandswan(1) = nobup
    elseif ( hasnobdown ) then
      numbandswan(1) = nobdown
    else 
      numbandswan(1) = noccupied(1)
      if( IOnode )          &
 &      write(6,'(/,a)')    &
 &        'number_bands_wannier: Automatic number of bands for wannierization'
    endif
    if ( numbandswan(1) .gt. no_u )                                   &
 &    call die( 'More bands required than available number of states. &
 &               Reduce the number of bands.')
  else ! Spin polarized calculation 
    if ( hasnob ) then
      if ( .not. ( hasnobup .or. hasnobdown )) then
        numbandswan(1) = nob
        numbandswan(2) = nob
      elseif ( .not. ( hasnobup .and. hasnobdown )) then 
        if( IOnode )          &
 &        write(6,'(a)')      &
 &        'number_bands_wannier: inconsistent output band specification'
          call die('Confusion in fdf')
      else 
        if( IOnode )          &
 &        write(6,'(a)')      &
 &        'number_bands_wannier: ignoring Siesta2Wannier90.NumberOfBands'
        numbandswan(1) = nobup
        numbandswan(2) = nobdown
      endif
    elseif ( hasnobup .and. hasnobdown ) then
      numbandswan(1) = nobup
      numbandswan(2) = nobdown
    elseif ( hasnobup .or. hasnobdown ) then
      if( IOnode )          &
 &      write(6,'(a)')      &
 &      'number_bands_wannier: inconsistent output band specification'
        call die('Confusion in fdf')
    else 
      if( IOnode )          &
 &      write(6,'(/,a)')    &
 &        'number_bands_wannier: Automatic number of bands for wannierization'
    endif
    if ( numbandswan(1) .gt. no_u .or. numbandswan(2) .gt. no_u )        &
 &    call die( 'More bands required than available number of states.    &
 &               Reduce the number of bands.')
  endif

  if( IOnode ) then   
    write(6,'(a)')                                                     &
 &    'number_bands_wannier: Number of bands for wannierization   '
    write(6,'(a,2i5)')                                                 &
 &    'number_bands_wannier: before excluding bands             = ',   &
 &     numbandswan
  endif

end subroutine number_bands_wannier

!
!-----------------------------------------------------------------------
subroutine set_excluded_bands( ispin, numexcluded, excludedbands, numbands, &
 &                             isexcluded, isincluded, numincbands )
!
! In this subroutine we set up the bands that are excluded from wannierization
!
! Used module procedures
  use alloc,           only: re_alloc ! Re-allocation routine
  use alloc,           only: de_alloc ! De-allocation routine

! Used module parameters
  use atomlist,        only: no_u     ! Number of orbitals in unit cell
  use parallel,        only: IOnode   ! Input/output node

!! For debugging
  use parallel,        only: Node     ! Working node
  use sys,             only: die      ! Termination routine
!! End debugging

  implicit none

  integer,  intent(in)            :: ispin
  integer,  intent(in)            :: numexcluded 
                                      ! Number of bands to exclude from the 
                                      !   from the calculation of the overlap 
                                      !   and projection matrices.
                                      !   This variable is read from the .nnkp 
                                      !   file
  integer,  pointer   :: excludedbands(:)
                                      ! Bands to be excluded.
                                      !   This variable is read from the .nnkp 
                                      !   file
  integer,  intent(in)            :: numbands(2)
  logical,  pointer   :: isexcluded(:)  ! List of bands excluded for
                                                    !   Wannierization
  integer,  pointer   :: isincluded(:)  ! List of bands included for
                                                    !   Wannierization
  integer,  intent(out)           :: numincbands(2)

!
! Internal variables
!
  integer :: iterex                   ! Counter for the loop on excluded bands
  integer :: iband                    ! Index of the band that is excluded
  
  nullify( isexcluded )
  call re_alloc( isexcluded, 1, no_u, name='isexcluded', &
 &               routine='set_excluded_bands' )

! By default, the occupied bands for the corresponding spin, numbands(ispin),
! are not excluded from the calculation...
  isexcluded(1:numbands(ispin))      = .false.

! ... while the unoccupied bands (from numbands(ispin+1) to the total number
! of bands (no_u)) are excluded
  isexcluded(numbands(ispin)+1:no_u) = .true.

  do iterex = 1, numexcluded
    iband = excludedbands( iterex )
    if ( iband .lt. 1 ) then
      call die("Erroneous list of excluded bands.")
    elseif( iband .gt. numbands(ispin) ) then
      call die("Erroneous list of excluded bands or numbands too small")
    endif
!   Exclude the corresponding band from the computation
    isexcluded( iband ) = .true.
  enddo

  call de_alloc( excludedbands, name='excludedbands',  &
 &               routine='set_excluded_bands')

! The number of actually "included" bands is the computation of 
! the overlap and projection matrices is the total number of 
! occupied bands minus the excluded bands
! This number will be used to dimension Mmnkb and Amnmat matrices,
! and also corrsponds to the number 
! of wave functions in UNK00001.1
  numincbands(ispin) = numbands(ispin) - numexcluded

! Determine the bands explicitly included for wannierization
  nullify( isincluded )
  call re_alloc( isincluded, 1, numincbands(ispin), name='isincluded', &
 &               routine='set_excluded_bands' )

  iband = 1
  do iterex = 1, numbands(ispin)
    if ( isexcluded( iterex ) .eqv. .false. ) then   
      isincluded( iband ) = iterex
      iband = iband + 1
    endif
  enddo 

  if( IOnode ) then   
    write(6,'(/,a,2i5)') 'Number of bands for wannierization ' // &
                         'after excluding bands:',  numincbands
    write(6,'(a)') 'Bands to be wannierized: '
    write(6,'(16i5)') (isincluded( iband ), iband=1, numincbands(ispin))
  endif
 
end subroutine set_excluded_bands

endmodule m_digest_nnkp
